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Abstract 

We present a time-saving simulator within the framework of the density functional theory to 
calculate the transport properties of electrons through nanostructures suspended between semi- 
infinite electrodes. By introducing the Fourier transform and preconditioning conjugate- gradient 
algorithms into the simulator, a highly efficient performance can be achieved in determining scat- 
tering wave functions and electron-transport properties of nanostructures suspended between semi- 
infinite jellium electrodes. To demonstrate the performance of the present algorithms, we study the 
conductance of metallic nanowires and the origin of the oscillatory behavior in the conductance of 
an Ir nanowire. It is confirmed that the s-d z 2 channel of the Ir nanowire exhibits the transmission 
oscillation with a period of two-atom length, which is also dominant in the experimentally obtained 
conductance trace. 

PACS numbers: 
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I. INTRODUCTION 



The understanding and control of the transport properties of electrons through nanos- 
tructures are important subjects for the development of new electronic devices. In this 
decade, numerous experimental and theoretical investigations Jiave been performed on the 

. studies, 



4] has 



atomic geometries and transport properties of nanostructures In the theoretica _ 
the nonequilibrium Green's function method within the tight-binding approach 
been widely employed to describe the quantum transport in nanoscale devices. While this 
method is efficient for investigating the transport properties of large systems, the complete- 
ness of a basis set is always a concern and the precise description of tunneling transport is 
hard work. 

To obtain exact theoretical knowledge on electron-transport properties, Fujimoto and 
Hirose presented the overbridging boundary-matching (OBM) method (5), 6|, formulated by 



the real-space finite-difference (RSFD) approach [6Ml l| within the framework of the density 
functional theory [l2, k|. Since the system is divided into the equally spaced grid points 
in the RSFD approach and the wave functions and potentials are directly defined on the 
grid points, one can strictly treat systems with arbitrary boundary conditions and eliminate 
difficulties arising from the incompleteness of a basis set. Therefore, the OBM method 

n 

enables us to examine tunneling transport [14[ as well as ballistic transport 



15-2J in 



nanostructures suspended between semi-infinite electrodes with a high degree of accuracy. 
However, in the OBM method, a large computational cost is required to calculate the Green's 
function necessary for determining the scattering wave functions. Recently, to avoid this 
computational hard work, Kong, Tiago, and Chelikowsky made an improvement to the 
OBM method 30|, enabling the scattering wave functions to be obtained by solving a set 



of simultaneous equations without the need to calculate the Green's function. We call 
this method the improved OBM (IOBM) method. However, the IOBM method is still 
inconvenient when electrodes have a large cross-sectional area in the transition region or are 
made of multivalent materials. 

In this paper, we propose efficient algorithms to solve the simultaneous equations aris- 
ing in the IOBM procedure and demonstrate electron-transport calculations for nanoscale 
junctions. To exemplify the advantages of our algorithms, we apply them to investigate the 
electron-transport properties of a single-row Na nanowire and find that the algorithms can 
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give reasonable numerical solutions within a short CPU time. We also examine the electron- 
transport properties of single-row Ir and Au monoatomic nanowires and explore why the 
even-odd oscillation is dominant in the experimentally obtained conductance traces. The 
results indicate that the even-odd oscillation is insensitive to the structural deformation of 
the nanowire, while oscillations with a longer period than two-atom length are easily affected 
by structural deformation. This implies that only the even-odd oscillation can survive and 
the other oscillation patterns are cancelled out in experiments. 

The remainder of this paper is organized as follows: Section |TT] gives the details of the 
computational scheme used to develop an efficient IOBM simulator to determine the electron- 
transport properties of nanostructures. In Sec. IIIIl the performance of our method is demon- 
strated. Furthermore, in Sec. IIV[ we adopt the method to examine the transport properties 
of Ir and Au atomic nanowire models and discuss the origin of the oscillatory behavior of 
conductance. Finally, our works are summarized in Sec. |V] and mathematical details are 
described in Appendix. 



II. COMPUTATIONAL FORMALISM 



We treat a system including a nanostructure suspended between semi-infinite electrodes as 
shown in Fig.[TJ Here, the x (y) and z coordinates are chosen to be parallel and perpendicular 
to the electrode surface, respectively, and the case of incident electrons propagating from 
the left electrode is considered. The system is infinite in the z direction and periodic in the 
x and y directions. The transition region is divided by grid points with equal spacing of 
hp = L^/N^ in the conventional OBM and IOBM methods based on the RSFD approach Jfj], 
where and iV M are the length and the number of grid points in the \i direction (/i = x, 
y, and z) of the transition region, respectively. The second-order derivative of the wave 
function ip in the Kohn-Sham equation is described as 

ip(r) « ^2 \ c Xy ni)(xi + nh x , y h z k ) + Cy >n if)(xi, y$ + nh y , z k ) + c g>n ij)(xi, y s , z k + nh z )j, 

(1) 



dr 2 

n=-N 



where Nf represents the parameter determining the order of the finite-difference approxima- 
tion and the coefficients c^ n are described in Ref. 7. The Kohn-Sham equation is written 
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in a discretized matrix form as 



B\*(z k - X ) + [E- H(z k ; «//)] (**) - B z ^{z k+1 ) = 

(k = — oo, 1,0,1,..., oo) 



(2) 



where B z = with I being the iV^-dimensional unit matrix (N xy = N x xN y ), H(z k ; K//) 

denotes the iV^-dimensional block-tridiagonal matrix including the potential on the x-y 
plane at z = z k , and nn = (k x , k v ) is the lateral Bloch wave vector within the first Brillouin 
zone. Here, *&(z k ) is the set of values of the wave function on the x-y plane at z — z k , 
{ip(xi,yj, z k ) : i = 1, • • • ,N X , j = 1, • • • ,N y }. For simplicity, the central finite-difference 
formula (Nf = 1) is adopted, and the extension to the case of a higher-order finite-difference 
representation is straightforward. In the transition region (zq < z < zn z +i), we rewrite 
Eq. d2J) in a matrix representation as 



E Hrp 



*(zo) 




' B^{z^) ' 


*(zi) 


















_B z q>(z Nz+2 )_ 



(3) 



where E is the Kohn-Sham energy and Ht is the Hamiltonian of the truncated part of the 
system sandwiched between the planes at z = z and z = z Nz+ i. H T forms the N x X N y X 
(N z + 2)-dimensional block-tridiagonal matrix given by 



H. 



T 



H(z ; K/i) B z 

B\ H(z 1 ;k // ) B z 











B\ H(z Nz ; K// ) B z 

B l #0jv*+i;k//) 



(4) 



Here and hereafter, the dependence of the Hamiltonian Ht on the lateral Bloch wave vector 
K/i is omitted from the suffix of variables to avoid complication. 

n n 

In the original OBM procedure |5|, [6[, the scattering wave functions are evaluated by 
employing the Green's function Qt{E) which is defined as the inverse matrix of [E — H T }- 
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When the block-matrix element of Qt{E) is denoted as Qij, Eq. (j3J) is rewritten as 



*(*b) 




£#0,0*(*-l) 


+ £*0O,JV.+l*(ZiV,+2) 






Bt0 lt(1 tt(z_i) 


+ J B z ^l,JV z+ l^'(2;Af z +2) 


tf(ZJV.+l) 




Bt&r, + i,o*(z-i) 


+ -B z ^Ar z + i iA r z+1 \l/(2;Ar z+ 2 



(5) 



Note that it is only necessary to work out the Oth and (N z + l)th block-column elements 
of Q T (E), Gift, and Gi,N x +i (i = 0, 1, • • • , JV 2 + 1). Since their calculation requires a compu- 
tational load proportional to 0(N% y x iV 2 ) with an iterative method such as the conjugate 
gradient (CG) method, a large computational load is required to adopt the OBM method 
to large systems. 



On the other hand, in the IOBM method 



30j , the scattering wave functions of electrons 



propagating from the left electrode are evaluated by solving the following simultaneous 
equations for each incident wave function with an iterative method: 



E-3 



T 



H 
























(6) 



Equation is obtained from Eq. ([3]) using the scattering boundary condition, i.e., 

Nxy 

t>?(z k )+J2ri®7 f (zk) (*<0) 
i=i 

N X y J 



(7) 



V 1=1 



where r$ (ij) is the reflection (transmission) coefficient, and $[ (;?&) ($- ra (zfc)) is the gener- 
alized Bloch state in semi-infinite electrodes for the reflected (transmitted) electrons. Here, 



H 



m*o) 










Y7 R (z Nz+ i) 



(8) 
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Note that Y7 L and YT R are the retarded self-energy matrices in the left and right electrodes, 
respectively, and are evaluated by 

E r L (z ) = BiQ^(z^)[Q^(z )}-\ (9) 

and 

Z R (z Nz+1 ) = B z Q tra {z Nz+2 )\Q tr \z Nz+l )}-\ (10) 

where Q A (z k ) (A = ref and tra) are iV^-dimensional matrices consisting of {<&f(zi~)}, i.e., 
Q A (z k ) = [$i(zk), • ■ ■ , ®N xv ( z k)] ( see Sec. 9.3 in Ref. 6). The generalized Bloch states 
inside the semi-infinite electrodes are constituted by propagating Bloch waves with real wave 
vectors and evanescent waves with complex ones. Because the evanescent waves behave as 
exponential functions and decay during their transmission from deep inside the left electrode, 
only right-propagating Bloch waves must be taken into account as incident waves The 
scattering wave functions for electrons propagating from the right electrode are also described 
in a similar manner. 

Since the retarded self-energy matrix is an N xy - dimensional matrix, the maximum order 
of the computational cost of solving Eq. (jHJ) is 0(N in x N% y ) owing to the multiplications 
of Y> r L (z ) x ty(zo) and T> t r (zn z +i) x ^/(zn z +i), where N in is the number of incident waves. 
Although [E — Ht — H] is a non-Hermitian matrix, we can save a reasonable amount of 
computational time compared with that required for the original OBM scheme when N in is 
much smaller than N z . On the other hand, as N in increases, this procedure for calculating 
the scattering wave functions consumes a larger CPU time than the original OBM scheme, 
because the convergence of the CG method for a non-Hermitian matrix is slow. Therefore, 
this procedure is not suitable when the system includes electrodes having a large cross- 
sectional area or consisting of multivalent materials, since N in is proportional to the lengths 
of sides of the supercell and the number of valence electrons within the electrodes. 

We now present a novel procedure for efficiently solving Eq. ([6]) in the case of jellium elec- 
trodes. The jellium-electrode approximation has been successfully app lied to the interpreta- 
tion of electron-transport properties with less computational load [15|, ll7|, Il9|, |2(J, 123J, |26|, 129 ]. 
The self-energy matrices of jellium electrodes are independent of Zk, and those in the left 
and right electrodes are the same, i.e., T, r L (z k ) = E^(z fc ) = Y7 . Thus, the component of the 
retarded self-energy matrix Y7 for the grid point (v/^V/,t') on the x-y plane in the case of 



Nf = 1 (the central finite-difference case) is analytically given by 



Sr 
e 01 



]— ex P + K //) • { r //,i - r //,£')] ■ exp(ik Zj „h z ), (11) 



(2ix 2tt \ 

where G// <v = (G Ux ,G Uy ) represents the lateral reciprocal lattice vector I — u x , —v y ) with 

\L X L y J 

v x and v y being integers, and k ZjV is the z component of the wave vector in the RSFD scheme 
defined as 



h 

"'Z.V 



COS ( — Xv) 

— for propagating waves {\Xv\ < 1) 

i cosh _1 (— x» 



(12) 



h. 



for evanescent waves (\xv\ > 1) 



with 



h 2 h 2 
X v = h 2 z E--^(l- cos(G Ux + Kx )h x ) - (1 - cos(G^ + K y )h y ) - 1. (13) 

a; y 

In Eq. (112j) . the branch of the function cos^ 1 ^) (cosh -1 (z)) is so chosen that k ZjU {—ik ZyU ) 
is a positive value. The derivation of Eqs. ( TTT]) -(fT3"l) is given in Appendix. 

The product of the matrix S r and the vector ty(zk), Y^jiipty/^r, z^), is in a convolution 

v 

form of the 2-dimensional Fourier transform, and it can be most easily calculated in the 
momentum space. Each matrix element of the Fourier transformed retarded self-energy 
matrix, Yt r uv ,, is given by 

Ky = JF- Yl Yl exp H(%* + ' V/A ex P H G //y + K ll) ■ Vi 'A ^le- ( 14 ) 

V/J V/,e' 

By applying the orthogonality of the plane waves, one obtains the following expression of 
the diagonalized self-energy matrix: 

Ely = -7^2$vu>exp{ik z yh z ). (15) 

The columnar vector ^(z^) consists of N xy values of the Fourier transformed wave function 
of ty(zk) as 

i'u(Zk) = — !- Y eXP M**//.* + K ll) ' V/A ^(Vf/^k)- ( 16 ) 

Because the number of grid points riii corresponds to that of the reciprocal lattice vectors 
Gji v and the off-diagonal elements of the Fourier transformed retarded self-energy matrix 



S r are zero according to Eq. (fT5|) . calculating the product of S r and ^(zk) requires a com- 
putational load of 0(N in x N xy ). The Fourier transform of ^{z^) and the inverse Fourier 
transform of S' r x ^(zk) require 0(Ni n x N xy logN xy ) operations when the fast Fourier 
transform (FFT) algorithm is employed. Thus, the overall scaling of the critical part of the 
calculations is improved from 0(N in x N% ) to 0(N in x N xy \ogN xy ). 

Next, we introduce a preconditioning CG (PCG) method into the solver to accelerate 
convergence since the total CPU time is proportional to the number of iterations. If the 
preconditioner P is similar to \E — — H\ , the spectral property of the matrix Fx [£- 
H T - H] can contribute to rapid convergence. The incomplete Cholesky preconditioners 
are effective at accelerating convergence and are commonly employed, but they require 
more computations per iteration and larger amounts of memory in general because the 
preconditioners are not sparse. One might consider taking an easily calculated and/or easily 
stored matrix P as an approximation to \E— Ht— H\ 1 . Although the Jacobi preconditioner, 
in which the preconditioner P is chosen to be the diagonal matrix in terms of the reciprocal 
of the diagonal elements of [E — Ht — H] , is one of the simplest and most useful forms, 
more sophisticated preconditioning is required for faster convergence. Here, we propose a 
method for improving P using the Green's function of the Laplacian operator. We postulate 
that the Laplacian of the kinetic operator in [£? — Ht — H\ is dominant and use the inverse 
matrix of the discretized Laplacian in the RSFD approach. Since the Green's function of 
the Laplacian operator is represented by -, the inverse matrix of the discretized 

Laplacian is approximated as - : , where raj) is the position of the i(j)th grid point. 

I r i ~ r j I 

To avoid the numerical difficulties due to this matrix not being sparse and diverging at 
Vi = rj, we employ the following truncated matrix as the preconditioner: 

Ay(r) = C o7 (r) + £ C l7 (r + r'), (17) 

r'=(±h x ,±h y ,±h z ) 

where 7 is the residual vector of the CG iteration, and Cq and C\ are coefficients used in 
the preconditioning, which are discussed later. We confirmed that this preconditioning can 
successfully reduce the CG iteration count required to solve the Poisson equation jf]]. 



S 



III. PERFORMANCE TEST 



In order to demonstrate the performance of the IOBM method incorporated with the FFT 
and PCG algorithms, the electron-transport properties of single-row Na nanowire models 
suspended between semi- infinite electrodes are examined (see Fig. [T]). The number of atoms 
consisting the nanowire, N at0 m, is varied between 3 and 8, and the nanowires are directly 
attached to structureless jellium electrodes. The interatomic distance is d (= \/3ao/2), where 
ao (=8.11 a.u.) is the lattice constant of Na bulk. The distance between the edge atom of the 
nanowire and the surface of the jellium electrode is A/2ao/4, so as to correspond to a [110] 
Na strand. The Wigner-Seitz radius of the jellium electrodes is taken to be r s = 3.99 a.u. 



The exchange-correlation effects are treated by the local density approximation 



31] and 



the interaction between electrons and atomic cores is described by the norm-conserving 



pseudopotentials of Troullier and Martins 32|, |33j . To determine the Kohn-Sham effective 



potential, a conventional supercell is employed under a periodic boundary condition in all 
directions, which is represented by a rectangle denoted by dashed lines in Fig. [TJ the lengths 
of the supercell are L x ^ y \ = 25.66 a.u. in the x(y) direction and L z = (N atam — 1)-g?+42.12 a.u. 
in the z direction. The numbers of grid points in the x(y) and z directions are set to be 
N x (y) = 36 and N z = 60 + 10 x (N atom — 1), respectively, and the number of incident 
waves, Ni n , is 13 at the Fermi level. In this case, the ratio of the maximum order for 
the calculation of E r x \P(zo(jv z +i)) with the FFT algorithm to that without adopting the 
Fourier transform is Nm * Nx y l °!= Nx y _ 0.0024. The conductance of the nanowire system 

■LV in J.\ X y 

at the limits of zero temperature and zero bias is determined by the Landauer-Biittiker 
formula j^J]. The numerical examination is carried out on a workstation with a 3.25 GHz 
Intel® Xeon® processor using four different solvers: (a) the original OBM method, (b) the 
IOBM method, (c) the IOBM method with the FFT algorithm, and (d) the IOBM method 
with the FFT and PCG algorithms. In solver (d), we adopt Cq = 1 and C\ = exp(— a) with 
a being a positive real number. 

CPU time versus N atom for solvers (a), (b), and (c) is displayed in Fig. [2j Note that 
solver (c) results in faster convergence than the others. To investigate the effect of the 
parameter of the PCG algorithm, Fig. [3] shows CPU time versus N atom for solver (d) with 
a = 3.0, 2.4, and 1.8. The CPU time of solver (c) is also plotted in Fig. [3] for comparison. 
Table [J shows the average number of iterations required to obtain convergent solutions for 
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incident waves. The average numbers of iterations for solver (d) with a = 3.0, 2.4, and 1.8 
are approximately 1.3, 1.7, and 2.3 times smaller than that for solver (c), respectively. The 
combination of the FFT and PCG algorithms enables us to significantly reduce the CPU 
time required for the calculation of the electron-transport properties of nanostructures with 
moderate memory consumption. 

Figure H] shows the conductance of Na nanowires in a unit of Go (Go = 2e 2 /h, where 
e is the electron charge and h is Planck's constant). The conductance is ~ 1 Go and 
exhibits oscillatory behavior with a period of two-atom length, i.e., the well-known even- 
odd oscillation. In previous theoretical studies, it was demonstrated that the conductance 
of a Na nanowire is not significantly different from 1 Go and that it oscillates with respect 



to the nanowire length [15 
of previous studies. 



3,Q, 



20 



27 



281 ] . Our result is in good agreement with those 



IV. APPLICATIONS 



In 2003, Smit et al. 21( found, using mechanically controllable break junctions, that 



aehavior with a 
4J carried out a 



the conductance of Ir and Pt monoatomic nanowires manifests oscillatory 
period of two-atom length, similarly to Na and Au ones, de la Vega et al. 
tight-binding calculation on the electron-transport properties of Ir, Pt, and Au nanowires 
and reported that additional oscillation patterns with a longer period and larger amplitude 
than those obtained in Ref. 2 can be observed in the conductance traces of Ir and Pt 
nanowires. Although it is intuitively expected that patterns with a large amplitude are 
dominant in the conductance traces, no experimental evidence of such oscillation patterns 
has been measured up to now. Recently, one of the present authors (T. O.) has examined 
the transport properties of a Pt nanowire by first-principles calculation and claimed that the 
even-odd oscillation is due to the low sensitivity of the transmission oscillation of the s-d z i 
channel to the spatial deformation of the nanowire [291 ] . On the other hand, in the case of the 
Ir nanowire, the transmission of the s-d z 2 channel oscillates with a longer period than that of 
the Pt nanowire according to a tight-binding calculation [4|. It is of interest to explore using 
first-principles calculations whether the even-odd oscillation of the conductance trace of the 
Ir nanowire is caused by the s-d z 2 channel. In addition, since 5d transition metals have the 
multiple valence electrons, Ni n for Ir electrodes is larger than that for Na ones, which causes 
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a significant increase in computational cost. For this reason, we apply the method described 
in the previous section to examine the transport properties of the Ir nanowire. 

We first calculate the electronic structures of infinite straight Ir and Au wires of equal 
interatomic distance d (= aj/y/2), where af is the lattice constant of fee bulk (a/ = 7.25 
and 7.71 a.u. for Ir and Au, respectively). The grid spacing h is set to be \/3af/34 and a 



double-grid technique [6 





5 is 


6. 
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condition, and the size of the supercell is L x = L y = 46h and L z = d, where L x (L y ) and L z 
are the lengths of the supercell in the x(y) directions perpendicular to the wire and in the z 
direction parallel to the wire, respectively. The exchange-correlation effects are treated by 
the local density approximation 31] and the interaction between electrons and atomic cores 



is described by the norm-conserving pseudopotentials of Troullier and Martins [32|, [33|]. The 
integration over the Brillouin zone along the wire-axis direction is performed by the sampling 
of 80 equidistant Appoints. We verified that the increase in the numbers of grid points and 
/c-points did not affect our conclusion. Figure [5] shows the energy band structures of the 
infinite Ir and Au wires. The bands are labeled according to the atomic orbitals mainly 
constituting the bands. Only the upper s-d z 2 band crosses the Fermi level in the case of 
the Au wire, whereas the other d bands also cross the Fermi level in the case of the Ir wire. 
These results are consistent with experimental results, in which the maximum conductance 
of the Au nanowire is 1 Go while the conductance of the Ir nanowire exceeds 1 Go 21]. 
In addition, it is well known that the oscillatory behavior of the conductance is led by the 
quantum-mechanical wave character of the electrons, and that the period of the oscillation 
is given by ir/k z at the intersection between the band and the Fermi level, where k z is the 



component of the wave vector along the wire axis in the infinite atomic wire [28j|. Since k z 
for the upper s-d z 2 channel is n/2d for both the Au wire and the Ir one, the transmission of 
these channels is expected to exhibit even-odd oscillation. 

We next examine the electron-transport properties of Ir and Au nanowires suspended 
between semi-infinite electrodes to ensure that the even-odd oscillation is observed in the 
transmission of the upper s-d z 2 channel. Hereafter, the superscripts Ir and Au denote the 
parameters for Ir and Au nanowires, respectively. The interatomic distances are the same 
as those of the infinite wires. The distance between the edge atom of the nanowire and the 
surface of the jellium electrode is 2o//v6. The Wigner-Seitz radii of the jellium electrodes 
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are taken to be r 1 / = 1.36 a.u. and r^ u = 1.35 a.u. so that the elements of the electrodes 
correspond to those of the nanowires. The grid spacing of d/16 is employed and N atom is 
varied between 3 and 8. The lengths of the transition region in the x and y directions 
are Li r = L\j = 15.07 a.u. and L^ u = L^ u = 15.38 a.u., and that in the z direction is 

x y x y i 

L z = (N atom + 4) • d a.u. The numbers of grid points in the x(y) direction are set to be 
Nxly) = 46 and N£j£\ = 48, and that in the z direction is N z = 16 x (N atom + 4). The number 
of incident waves, iV in , is 37 for both models. In this case, the ratios of the maximum 
order for the calculation of S r x ^{zo(n z +i)) with the FFT algorithm to that without the 
Fourier transform are = 0.0016 and = 0.0015 for the Ir and 

Au nanowires, respectively. The other computational conditions are the same as those for 
the infinite wires. 

Figure [6] shows the conductance of the Ir and Au nanowires in the unit of Go- The 
conductance trace of the Au nanowire exhibits even-odd oscillation depending on N atom 
with an amplitude of ~ 0.02 Go, while no oscillatory behavior is observed in the conductance 
trace of the Ir nanowires. To obtain a deeper understanding of the oscillatory behavior, the 
evolution of the channel decomposition of the conductance trace is shown in Table HIl in 
which the quantum numbers of the eigenchannels correspond to those in the case of infinite 

n 

wires. The eigenchannels are computed by diagonalizing the Hermitian matrix T^T [171 ]. 
where T is the transmission matrix. Even-odd oscillation clearly emerges in the transmission 
of the upper s-d z i channel of the Ir nanowire, and the other d channels exhibit oscillation 
with a longer period and a larger amplitude of ~ 0.5 Go- According to a previous theoretical 
study on Pt nanowires 29( , the contributions of the other d channels become negligible upon 
being averaged over thousands of scans in experiments. Thus, we can conclude that the even- 



odd oscillation in the conductance trace of the Ir nanowire observed in the experiment 
is attributed to the transmission oscillation of the upper s-d z 2 channel. 



2l| 



V. CONCLUSION 



We have presented algorithms for efficiently solving a set of simultaneous equations de- 
scribed in the IOBM method proposed by Kong, Tiago, and Chelikowsky [30j. In the case of 
structureless jellium electrodes, the self-energy matrices are analytically given and are diag- 
onal in momentum space. These spectral properties of the self-energy matrices enable us to 
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reduce the computational cost required for the multiplication of the self-energy matrices and 
wave functions. Furthermore, a solver that successfully accelerated the CG algorithm has 
been developed by introducing a preconditioning technique. In this way, we developed an 
extremely efficient simulator. To demonstrate the performance of the presented algorithms, 
we applied them to calculate the transport properties of an Ir nanowire attached to semi- 
infinite electrodes. The upper s-d z 2 channel of the Ir nanowire mainly contributes to the 
electron transport and gives rise to the oscillation of the transmission with a period of two- 
atom length. Our accelerated algorithms are expected to lead to a greater understanding of 
the physics underlying electron transport through nanostructures using large computational 
models. 
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Appendix 



Retarded self-energy matrix of a semi-infinite jellium electrode 

As shown in Eqs. and f lTU]) . the retarded self-energy matrices Y7 L and are obtained 
by using the generalized Bloch states Q re ^(zk) and Q tra (zk)- On the other hand, the self- 
energy matrices are also derived from the retarded Green's function Q^ L R ^ in the semi-infinite 
electrodes (see p. 160 of Ref. 6) as 

££0o) = B\QI{z^z^E,ku)B z , (Al) 

and 

Zr( z n z+ i) = B z g r R {z Nz+2 ,z N ^2,E,K//)Bl. (A.2) 

In this appendix, we introduce the analytical derivation of the retarded Green's function 
in the left-hand semi-infinite jellium electrode. That in the right-hand one is derived in 
the same manner. For simplicity, the lateral Bloch vector Kn is taken to be zero. We 
deal with the model of a free electron in which the discretized space is semi-infinite in 
the z direction and periodic in the x and y directions (see Fig. [7]). The grid points are 
denoted by r t = {x ix ,yi u , z h ) = (l x h x ,l y h y ,l z h z ), where Z M = -iV^/2, -1,0, 1, A^/2-l, 
and and are the grid spacing and the total number of the grid points in the /i 
direction (fi = x,y), respectively. We first assume a finite system in the z direction, i.e., 
l z = —N z — 1, — N z , • • • , —1, and impose the zero boundary condition on the wave function, 
4>m(z-N z -i) = <Pm{zo) — 0, and then a semi-infinite system will be represented by taking 
the limit of N z — > oo. Adopting the central finite-difference formula for the second-order 

derivative, we write the Schrodinger equation for a free electron according to Eq. ([[]) as 
j i 

- g { c ^^(pm(xi x + nh x , y iy , z h ) + Cy tn (f> m (xi x ,y ly + nh y , z h ) 

n=—l 

+c z , n <frm(xi x ,y ly ,z lz +nh z )\ = E m (f) m (xi x ,yi y ,z lz ) 

(A.3) 

with c^o = —2/h and c M) ±i = (fi = x,y, z). In a finite system, the Hamiltonian has 
discrete eigenvalues and therefore the Green's function 

Q(r u r v -Z) = > (A.4) 
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exhibits only simple poles at the positions of the eigenvalues in the complex Z plane, where 
{E m } is the set of the eigenvalues and {4>m(r)} is the complete orthonormal set of the 
eigenf unctions of the Hamiltonian. In this model, the eigenfunctions and eigenvalues are 
analytically described as 



2 



m(xi x ,yi y ,zi z ) = J N — + ^ exp{i(G Vx x lx +G Uy yi y )j sin $ Vz z lx (A.5) 



and E m = ^-(1 - cos G Vx h x ) + — (1 - cos G Uy h y ) + — (1 - cos (3 Uz h z ), (A. 6) 



71 

respectively, where j3 Uz = — ^ v z and N xy — N x x N y . Thus, the Green's function of 
this finite system is given by 



n x i N y 

2 x 2 



2 , 

Q L {r u r v -Z) = N ( Nz + ^ X, exp{iG^(a; fc - ^)} exp{iG^(^ -y^)} 



— a" "w= — r 



sin/^z;, sin/3^/ 

X ^ Z - 1 ' j 

v x =\ 

We go on to the derivation of the Green's function for a free electron in the semi-infinite 
discretized space. By carrying out a limiting procedure N z — > oo while keeping h z constant 
in ( 1A.7I) . we have 



N X 1 Ny 

h 2 2 2 

g L (r h r V ]Z) = ^ ^ expjiG^^ - x z /J} expjiG^^ -?/«;)} 

XJ/ ,. _ N X Ny 

v x - — 2" v y =—f 

x ; exp{^(/ 2 -/;)}-exp{^(/ z + ^)} ^ 

-7T /^Z - {^1(1 - cosG Vx h x ) + jjf (1 - cosG Vy h y ) + (1 - cos#)} 



(A.8) 



To evaluate the integral in Eq. (1A.8j) , we transform it to an integral over the complex variable 
oj (= e ) along the unit circle as 

N x -, JV;; -i 
^2 2 2 

Q L {r h r v \Z) = — Yl e ^{iG Vx {x lx - x Vx )} &vp{iG Vy (y ly -y Vy )} 
mi\ xy _ 

Vx- 2" Vy=--f- 

J\u>\=i - wi(z)J (w - o; 2 (z) ) 
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where 



LU 2 (Z) = - X u{Z) - y/ X u(Z) 2 -l 

X„(Z) = h\Z - ^ (1 - cosGM - ^§(1 - cosG^/g - 1. 



(A.10) 



It follows that bj\bJi = 1. Only when the pole u = u\ or W2 in ( 1A.9I) exists inside the unit 
circle in the u plane, does it contribute to the integral. In a similar manner as shown in 
p. 149-151 of Ref. 6, the retarded Green's function is derived as 



Ql{r h r v -E) 



h 



— i Sk-i 

2 2 2 



iN, 



S exp{iG Vx {xi x - xv x )} exp{iG Vy {yi y - y Vy )} 



Xy „ - N x 

fx- — Uy= f 



exp{ik Zt „\zi z - zu\} - exp{ik ZtU \zi z + z v \} 



(All) 



sin k z>v h z 

with k z>u defined by Eq. (fl2l) . Finally, the retarded Green's function at the surface of the 
left-hand semi-infinite electrode is represented by 



Gl i x i* , x i' x ,yi y ,yv v ,z-uz-x,] E) 



2h 



N X 1 Ny 

2 2 2 



A 7 



J2 exp^'G^fa. -x Vx )} 



x y ,. _ n x 



2 "v— — 



x exp{iG Uy (y ty - y Vy )} exp(ik z , v h z ) . (A.12) 

The above-mentioned derivation of the Green's function is straightforwardly extended to the 
case of a higher-order finite-difference approximation. 
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FIG. 1: Schematic view of ballistic transport system. The transition region is suspended between 
the left and right semi-infinite electrodes. In the x and y directions, periodic boundary conditions 
are imposed. The dashed lines represent the boundaries of the supercell employed to determine 
the Kohn-Sham effective potential and the transition region employed to calculate the transport 
properties. The shaded area represents jellium electrodes. 




FIG. 2: CPU time required to calculate the conductances of single-row Na nanowires as a function 
of the number of atoms constituting the nanowires, N a t m- (a) original OBM method, (b) IOBM 
method, and (c) IOBM method with the FFT algorithm. 
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FIG. 3: CPU time required to calculate the conductances of single-row Na nanowires as a function 
of the number of atoms constituting the nanowires, N a t om , employing the IOBM method with the 
FFT and PCG algorithms, (c) and (d) correspond to the solvers described in the text. 
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FIG. 4: Conductance of single-row Na nanowires as a function of the number of atoms constituting 
the nanowires, N atom . 
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FIG. 5: Energy band structures of infinite (a) Ir and (b) Au wires. The zero of energy is chosen 
to be the Fermi level. 
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FIG. 6: Conductances of (a) Ir and (b) Au nanowires as a function of the number of atoms 
constituting the nanowires, N at0 m- 
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FIG. 7: Schematic view of discretized 3-dimensional semi-infinite systems. 
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TABLE I: Average CG iteration counts required to calculate scattering wave functions for inci- 
dent waves in single-row sodium nanowires employing the IOBM method with the FFT and PCG 
algorithms. 



Aatom 


without PCG 


a = 3.0 


a = 2.4 


a = 1.8 


3 


13936.46 


10864.38 


8241.92 


6107.30 


4 


13852.77 


10738.76 


8184.76 


6030.84 


5 


14688.00 


11518.92 


8918.15 


6620.00 


6 


14058.38 


10995.69 


8336.00 


6326.38 


7 


14639.15 


12036.07 


9111.92 


6750.23 


8 


14818.23 


11597.76 


8794.30 


6613.76 



TABLE II: Channel transmissions of Ir nanowires. 



n 


upper s-d z 2 


lower s-d z 2 




dx 2 — y 2 (xy) 


3 


0.854 


0.016 


0.763 


0.016 


4 


0.973 


0.011 


0.935 


0.497 


5 


0.727 


0.018 


0.510 


0.019 


6 


1.000 


0.020 


0.656 


0.010 


7 


0.838 


0.056 


0.979 


0.048 


8 


0.975 


0.116 


0.569 


0.206 
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